#ifndef SPLINE_HPP
#define SPLINE_HPP

#include <stdexcept>
#include <cmath>
#include <vector>
#include <iostream>
#include <cblas.h> 
#include <lapacke.h> 
#include "Function.hpp"

class Spline {
protected:
    std::vector<double> node;  // 样条节点
public:
    Spline(const std::vector<double>& node) : node(node) {}

    virtual double getSplineValue(double x) const = 0;  

    virtual ~Spline() = default;  // 默认析构函数
};

class LinearPPForm : public Spline {
private:
    std::vector<double> y_vals;   // Store y values
    std::vector<double> y_prime;  // Store first derivatives (slopes)

    // Set first derivatives based on linear interpolation
    void setFirstDerivatives(const std::vector<double>& x_vals, const std::vector<double>& y_vals) {
        int n = x_vals.size();
        y_prime.resize(n - 1);

        // Linear interpolation: the slope between each pair of consecutive points
        for (int i = 0; i < n - 1; ++i) {
            y_prime[i] = (y_vals[i + 1] - y_vals[i]) / (x_vals[i + 1] - x_vals[i]);
        }
    }

public:
    // Constructor for LinearPPForm
    LinearPPForm(const std::vector<double>& x_vals, const std::vector<double>& y_vals)
        : Spline(x_vals), y_vals(y_vals) {
        // Compute the first derivatives
        setFirstDerivatives(x_vals, y_vals);
    }

    // Get the spline value at a given x using linear interpolation
    double getSplineValue(double x) const override {
        int n = node.size();
        int i = 0;

        // Find the interval that contains x
        for (i = 0; i < n - 1; ++i) {
            if (x >= node[i] && x <= node[i + 1]) break;
        }

        // Linear interpolation formula
        double h = node[i + 1] - node[i]; // Interval width
        double slope = y_prime[i];  // Slope (first derivative) between node[i] and node[i+1]
        double intercept = y_vals[i] - slope * node[i];  // Line equation: y = slope * x + intercept

        return slope * x + intercept;  // Linear interpolation
    }
};

class PPForm : public Spline {
private:
    std::vector<double> y_vals;   // 存储y值
    std::vector<double> y_prime;  // 存储一阶导数值
    std::vector<double> diff_div; // 存储差商

    // 将一阶导数数组拷贝到y_prime中
    void setFirstDerivatives(const double deriv[], int n) {
        y_prime = std::vector<double>(deriv, deriv + n);
    }

    // 计算自然边界条件
    void calcNaturalBoundary(const std::vector<double>& x_vals, const std::vector<double>& y_vals) {
        int n = x_vals.size();
        std::vector<double> lower_coeffs(n-1, 0.0), upper_coeffs(n-1, 0.0), diag_coeffs(n, 2.0), rhs(n, 0.0);

        // 设置自然边界：s''(a)=0, s''(b)=0
        lower_coeffs[n-2] = 1.0;
        for (int i = 0; i < n-2; ++i) {
            lower_coeffs[i] = (x_vals[i+3] - x_vals[i+2]) / (x_vals[i+3] - x_vals[i+1]);
        }

        upper_coeffs[0] = 1.0;
        for (int i = 1; i < n-1; ++i) {
            upper_coeffs[i] = (x_vals[i+1] - x_vals[i]) / (x_vals[i+2] - x_vals[i]);
        }

        rhs[0] = 3 * diff_div[0];  
        rhs[n-1] = 3 * diff_div[n-2]; 
        for (int i = 1; i < n-1; ++i) {
            rhs[i] = 3 * upper_coeffs[i] * diff_div[i] + 3 * lower_coeffs[i-1] * diff_div[i-1];
        }

        // 求解三对角方程得到一阶导数
        y_prime = solveTridiagonal(lower_coeffs, diag_coeffs, upper_coeffs, rhs);
    }

    // 计算周期边界条件
    void calcPeriodicBoundary(const std::vector<double>& x_vals, const std::vector<double>& y_vals) {
        int n = x_vals.size();
        std::vector<double> matrix_A(n * n, 0.0), rhs_vals(n, 0.0);

        std::vector<double> lower_coeffs(n-1, 0.0), upper_coeffs(n-1, 0.0), diag_coeffs(n, 0.0);

        // 构造周期性三对角矩阵
        diag_coeffs[0] = 2.0 / (x_vals[1] - x_vals[0]);
        diag_coeffs[n-1] = 2.0 / (x_vals[n-1] - x_vals[n-2]);
        for (int i = 1; i < n-1; ++i) {
            diag_coeffs[i] = 2.0;
            lower_coeffs[i-1] = (x_vals[i] - x_vals[i-1]) / (x_vals[i+1] - x_vals[i-1]);
            upper_coeffs[i] = (x_vals[i+1] - x_vals[i]) / (x_vals[i+1] - x_vals[i-1]);
        }

        // 填充常数项
        rhs_vals[0] = 3 * diff_div[0] / (x_vals[1] - x_vals[0]) + 3 * diff_div[n-2] / (x_vals[n-1] - x_vals[n-2]);
        for (int i = 1; i < n-1; ++i) {
            rhs_vals[i] = 3 * (x_vals[i] - x_vals[i-1]) / (x_vals[i+1] - x_vals[i-1]) * diff_div[i]
                        + 3 * (x_vals[i+1] - x_vals[i]) / (x_vals[i+1] - x_vals[i-1]) * diff_div[i-1];
        }

        rhs_vals[0] += 3 * diff_div[n-1] / (x_vals[1] - x_vals[0]);
        rhs_vals[n-1] += 3 * diff_div[0] / (x_vals[n-1] - x_vals[n-2]);

        // 求解三对角线性方程组
        y_prime = solveTridiagonal(lower_coeffs, diag_coeffs, upper_coeffs, rhs_vals);
    }

    // 完全边界条件
    void calcCompleteBoundary(const std::vector<double>& x_vals, const std::vector<double>& y_vals, const double first_der_a, const double first_der_b) {
        int n = x_vals.size();
        std::vector<double> lower_coeffs(n-3, 0.0), upper_coeffs(n-3, 0.0), diag_coeffs(n-2, 2.0), rhs_vals(n-2, 0.0);

        // 设置三对角矩阵
        for (int i = 0; i < n-3; ++i) {
            lower_coeffs[i] = (x_vals[i+3] - x_vals[i+2]) / (x_vals[i+3] - x_vals[i+1]);
        }
        for (int i = 0; i < n-3; ++i) {
            upper_coeffs[i] = (x_vals[i+1] - x_vals[i]) / (x_vals[i+2] - x_vals[i]);
        }

        rhs_vals[0] = 3 * upper_coeffs[0] * diff_div[1] + 3 * (x_vals[2] - x_vals[1]) / (x_vals[2] - x_vals[0]) * diff_div[0] - (x_vals[2] - x_vals[1]) / (x_vals[2] - x_vals[0]) * first_der_a;
        rhs_vals[n-3] = 3 * (x_vals[n-2] - x_vals[n-3]) / (x_vals[n-1] - x_vals[n-3]) * diff_div[n-2] + 3 * lower_coeffs[n-4] * diff_div[n-3] - (x_vals[n-2] - x_vals[n-3]) / (x_vals[n-1] - x_vals[n-3]) * first_der_b;

        for (int i = 1; i < n-3; ++i) {
            rhs_vals[i] = 3 * upper_coeffs[i] * diff_div[i+1] + 3 * lower_coeffs[i-1] * diff_div[i];
        }

        // 求解三对角线性方程组
        y_prime = solveTridiagonal(lower_coeffs, diag_coeffs, upper_coeffs, rhs_vals);
        y_prime.insert(y_prime.begin(), first_der_a);
        y_prime.push_back(first_der_b);
    }

    // 计算导数值
    void calcFirstDerivatives(const std::vector<double>& x_vals, const std::vector<double>& y_vals, const std::string& boundary_type, const double first_der_a, const double first_der_b) {
        if (boundary_type == "natural") {
            calcNaturalBoundary(x_vals, y_vals);
        } else if (boundary_type == "periodic") {
            calcPeriodicBoundary(x_vals, y_vals);
        } else if (boundary_type == "complete") {
            calcCompleteBoundary(x_vals, y_vals, first_der_a, first_der_b);
        }
    }

    void calcFirstDerivatives(const std::vector<double>& x_vals, const std::vector<double>& y_vals, const std::string& boundary_type) {
        if (boundary_type == "natural") {
            calcNaturalBoundary(x_vals, y_vals);
        } else if (boundary_type == "periodic") {
            calcPeriodicBoundary(x_vals, y_vals);
        }
    }

public:
    PPForm(const std::vector<double>& x_vals, const std::vector<double>& y_vals, const std::string& boundary_type, const double first_der_a, const double first_der_b)
        : Spline(x_vals), y_vals(y_vals) {
        int n = x_vals.size();
        diff_div.resize(n - 1);
        
        // 计算差商
        for (int i = 0; i < n-1; ++i) {
            diff_div[i] = (y_vals[i + 1] - y_vals[i]) / (x_vals[i + 1] - x_vals[i]);
        }

        // 根据边界类型计算一阶导数
        calcFirstDerivatives(x_vals, y_vals, boundary_type, first_der_a, first_der_b);
    }

    PPForm(const std::vector<double>& x_vals, const std::vector<double>& y_vals, const std::string& boundary_type)
        : Spline(x_vals), y_vals(y_vals) {
        int n = x_vals.size();
        diff_div.resize(n - 1);
        
        // 计算差商
        for (int i = 0; i < n-1; ++i) {
            diff_div[i] = (y_vals[i + 1] - y_vals[i]) / (x_vals[i + 1] - x_vals[i]);
        }

        // 根据边界类型计算一阶导数
        calcFirstDerivatives(x_vals, y_vals, boundary_type);
    }

    // 根据计算的导数值插值
    double getSplineValue(double x) const override {
        int n = node.size();
        int i = 0;
        
        // 查找包含x的区间
        for (i = 0; i < n - 1; ++i) {
            if (x >= node[i] && x <= node[i + 1]) break;
        }
        
        // 计算插值
        double h = node[i + 1] - node[i];

        // 使用一阶导数进行插值
        double a = y_vals[i];
        double b = y_prime[i];
        double c = (3 * diff_div[i] - 2 * y_prime[i] - y_prime[i + 1]) / h;
        double d = (y_prime[i] + y_prime[i + 1] - 2 * diff_div[i]) / (h * h);
        
        return a + b * (x - node[i]) + c * std::pow(x - node[i], 2) + d * std::pow(x - node[i], 3);
    }
};

class LinearBSpline : public Spline {
private:
    int spline_degree; // 样条阶数，固定为 1
    std::vector<double> knots;  // 节点位置
    std::vector<double> spline_values;  // 样条函数值
    std::vector<double> coefficients;  // 样条系数

    // 计算线性 B 样条基函数
    double B(int i, int degree, double t) const {
        if (degree == 0) {
            return (knots[i] <= t && t < knots[i + 1]) ? 1.0 : 0.0;
        }
        double left = (t - knots[i]) / (knots[i + degree] - knots[i]);
        double right = (knots[i + degree + 1] - t) / (knots[i + degree + 1] - knots[i + 1]);
        return left * B(i, degree - 1, t) + right * B(i + 1, degree - 1, t);
    }

    // 设置自然边界条件
    void applyNaturalBoundary() {
        int num_internal_knots = knots.size() - 2; // 内部节点数
        coefficients.resize(num_internal_knots, 0.0);
        for (int i = 0; i < num_internal_knots; ++i) {
            coefficients[i] = spline_values[i]; // 对于线性 B 样条，直接等于采样值
        }
    }

public:
    LinearBSpline(const std::vector<double>& x, const std::vector<double>& y, const std::string& boundary_type, double fa = 0, double fb = 0)
        : Spline(x), spline_values(y), knots(x), spline_degree(1) {
        // 插值节点的拓展，适应 B 样条计算
        int N = x.size();
        double length = x[N - 1] - x[0];
        knots.insert(knots.begin(), x[0] - length / (N - 1));
        knots.push_back(x[N - 1] + length / (N - 1));

        if (boundary_type == "natural") {
            applyNaturalBoundary();
        }
    }

    LinearBSpline(const std::vector<double>& x, const std::vector<double>& y)
        : LinearBSpline(x, y, "natural") {} // 默认使用自然边界条件

    double getSplineValue(double x) const override {
    // 检查 x 是否在有效范围内
    if (x < knots[spline_degree] || x > knots[knots.size() - spline_degree - 1]) {
        throw std::out_of_range("x is out of range.");
    }

    // 计算结果
    double result = 0.0;

    // 遍历有效的 B 样条基函数
    for (size_t i = 0; i < coefficients.size(); ++i) {
        result += coefficients[i] * B(i, spline_degree, x);
    }
    return result;
    }
};

class BSpline : public Spline {
private:
    int spline_degree;
    std::vector<double> knots;  // 节点位置
    std::vector<double> spline_values;  // 样条函数值
    std::vector<double> coefficients;  // 样条系数

    // 计算自然边界条件
    void applyNaturalBoundary() {
        int total_points = knots.size();
        int num_internal_knots = total_points - 1 - spline_degree;
        double matrix[num_internal_knots * num_internal_knots];
        double rhs[num_internal_knots];
        std::fill(matrix, matrix + num_internal_knots * num_internal_knots, 0.0);
        std::fill(rhs, rhs + num_internal_knots, 0.0);
        int pivots[num_internal_knots];
        int status;

        // 自然边界条件：二阶导数在端点处为零
        for (int i = 0; i < num_internal_knots; ++i) {
            matrix[i] = B_derivative(i, spline_degree, knots[spline_degree], 2);
        }

        // 填充内部点的矩阵
        for (int i = 1; i < num_internal_knots - 1; ++i) {
            for (int j = 0; j < num_internal_knots; ++j) {
                matrix[i * num_internal_knots + j] = B(j, spline_degree, knots[spline_degree + i - 1]);
            }
            rhs[i] = spline_values[i - 1];
        }

        // 设置最后一行，满足自然边界条件
        for (int i = 0; i < num_internal_knots; ++i) {
            matrix[num_internal_knots * num_internal_knots - num_internal_knots + i] = B_derivative(i, spline_degree, knots[num_internal_knots], 2);
        }

        // 使用 LAPACKE 解线性方程组 Ax = b
        status = LAPACKE_dgesv(LAPACK_ROW_MAJOR, num_internal_knots, 1, matrix, num_internal_knots, pivots, rhs, 1);

        if (status == 0) {
            coefficients = std::vector<double>(rhs, rhs + num_internal_knots);
        } else {
            std::cerr << "Error solving the system: " << status << std::endl;
        }
    }

    // 计算周期边界条件
    void applyPeriodicBoundary() {
        int total_points = knots.size();
        int num_internal_knots = total_points - 1 - spline_degree;
        double matrix[num_internal_knots * num_internal_knots];
        double rhs[num_internal_knots];
        int pivots[num_internal_knots];
        int status;

        std::fill(matrix, matrix + num_internal_knots * num_internal_knots, 0.0);
        std::fill(rhs, rhs + num_internal_knots, 0.0);

        // 设置周期样条条件：s(a) = s(b), s'(a) = s'(b), s''(a) = s''(b)
        for (int i = 0; i < num_internal_knots; ++i) {
            matrix[i] = B_derivative(i, spline_degree, knots[spline_degree], 2) - B_derivative(i, spline_degree, knots[num_internal_knots], 2);
        }

        // 填充内部点的矩阵
        for (int i = 1; i < num_internal_knots - 1; ++i) {
            for (int j = 0; j < num_internal_knots; ++j) {
                matrix[i * num_internal_knots + j] = B(j, spline_degree, knots[spline_degree + i - 1]);
            }
            rhs[i] = spline_values[i - 1];
        }

        // 处理边界条件
        for (int i = 0; i < num_internal_knots; ++i) {
            matrix[num_internal_knots * num_internal_knots - num_internal_knots + i] = B_derivative(i, spline_degree, knots[spline_degree], 1) - B_derivative(i, spline_degree, knots[num_internal_knots], 1);
        }

        // 使用 LAPACKE 求解 Ax = b
        status = LAPACKE_dgesv(LAPACK_ROW_MAJOR, num_internal_knots, 1, matrix, num_internal_knots, pivots, rhs, 1);

        if (status == 0) {
            coefficients = std::vector<double>(rhs, rhs + num_internal_knots);
        } else {
            std::cerr << "Error solving the system: " << status << std::endl;
        }
    }

    // 计算完整的边界条件
    void applyCompleteBoundary(const double first_derivative_start, const double first_derivative_end) {
        int total_points = knots.size();
        int num_internal_knots = total_points - 1 - spline_degree;
        double matrix[num_internal_knots * num_internal_knots];
        double rhs[num_internal_knots];
        std::fill(matrix, matrix + num_internal_knots * num_internal_knots, 0.0);
        std::fill(rhs, rhs + num_internal_knots, 0.0);
        int pivots[num_internal_knots];
        int status;

        // 完整的边界条件：s'(start) = first_derivative_start, s'(end) = first_derivative_end
        for (int i = 0; i < num_internal_knots; ++i) {
            matrix[i] = B_derivative(i, spline_degree, knots[spline_degree], 1);
        }

        // 填充矩阵 M
        for (int i = 1; i < num_internal_knots - 1; ++i) {
            for (int j = 0; j < num_internal_knots; ++j) {
                matrix[i * num_internal_knots + j] = B(j, spline_degree, knots[spline_degree + i - 1]);
            }
        }

        // 填充边界条件的矩阵
        for (int i = 0; i < num_internal_knots; ++i) {
            matrix[num_internal_knots * num_internal_knots - num_internal_knots + i] = B_derivative(i, spline_degree, knots[num_internal_knots], 1);
        }

        // 设置常数向量 b
        rhs[0] = first_derivative_start;
        for (int i = 1; i < num_internal_knots - 1; ++i) {
            rhs[i] = spline_values[i - 1];
        }
        rhs[num_internal_knots - 1] = first_derivative_end;

        // 使用 LAPACKE 求解 Ax = b
        status = LAPACKE_dgesv(LAPACK_ROW_MAJOR, num_internal_knots, 1, matrix, num_internal_knots, pivots, rhs, 1);

        if (status == 0) {
            coefficients = std::vector<double>(rhs, rhs + num_internal_knots);
        } else {
            std::cerr << "Error solving the system: " << status << std::endl;
        }
    }

    // 计算 B 样条基函数
    double B(int i, int degree, double t) const {
        if (degree == 0) {
            return (knots[i] <= t && t < knots[i + 1]) ? 1.0 : 0.0;
        }
        double left = (t - knots[i]) / (knots[i + degree] - knots[i]);
        double right = (knots[i + degree + 1] - t) / (knots[i + degree + 1] - knots[i + 1]);
        return left * B(i, degree - 1, t) + right * B(i + 1, degree - 1, t);
    }

    // 计算 B 样条导数
    double B_derivative(int i, int degree, double t, int order) const {
        double left = degree / (knots[i + degree] - knots[i]);
        double right = degree / (knots[i + degree + 1] - knots[i + 1]);
        if (order == 0) {
            return B(i, degree, t);
        }
        return left * B_derivative(i, degree - 1, t, order - 1) - right * B_derivative(i + 1, degree - 1, t, order - 1);
    }

public:
    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree, const std::string& boundary_type, const double fa, const double fb)
        : Spline(x), spline_values(y), spline_degree(degree), knots(x) {
        int N = x.size();
        double length = x[N - 1] - x[0];
        for (int i = 0; i < degree; ++i) {
            knots.insert(knots.begin(), x[0] - (i + 1) * length / (N - 1));
            knots.push_back(x[N - 1] + (i + 1) * length / (N - 1));
        }
        if (degree == 3) {
            if (boundary_type == "natural") {
                applyNaturalBoundary();
            } else if (boundary_type == "periodic") {
                applyPeriodicBoundary();
            } else if (boundary_type == "complete") {
                applyCompleteBoundary(fa, fb);
            }
        }
    }

    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree, const std::string& boundary_type)
        : Spline(x), spline_values(y), spline_degree(degree), knots(x) {
        int N = x.size();
        double length = x[N - 1] - x[0];
        for (int i = 0; i < degree; ++i) {
            knots.insert(knots.begin(), x[0] - (i + 1) * length / (N - 1));
            knots.push_back(x[N - 1] + (i + 1) * length / (N - 1));
        }
        if (boundary_type == "natural") {
            applyNaturalBoundary();
        } else if (boundary_type == "periodic") {
            applyPeriodicBoundary();
        }
    }

    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree)
        : Spline(x), spline_values(y), spline_degree(degree), knots(x) {
        int N = x.size();
        double length = x[N - 1] - x[0];
        for (int i = 0; i < degree; ++i) {
            knots.insert(knots.begin(), x[0] - (i + 1) * length / (N - 1));
            knots.push_back(x[N - 1] + (i + 1) * length / (N - 1));
        }
        applyNaturalBoundary(); // 默认自然边界
    }

    BSpline(const std::vector<double>& x, const std::vector<double>& y, int degree, const std::vector<double>& coefficients)
        : Spline(x), spline_values(y), spline_degree(degree), knots(x), coefficients(coefficients) {
        int N = x.size();
        double length = x[N - 1] - x[0];
        for (int i = 0; i < degree; ++i) {
            knots.insert(knots.begin(), x[0] - (i + 1) * length / (N - 1));
            knots.push_back(x[N - 1] + (i + 1) * length / (N - 1));
        }
    }

    double getSplineValue(double x) const override {
        int total_points = knots.size();
        int num_internal_knots = total_points - 1 - spline_degree;
        double result = 0.0;

        // 计算基函数并加权
        for (int i = 0; i < num_internal_knots; ++i) {
            result += coefficients[i] * B(i, spline_degree, x);
        }
        return result;
    }
};
#endif